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ABSTRACT 

As emphasized by previous studies, proper treatment of the density fluctuation on the fundamental 
scale of a cosmological simulation volume - the "DC mode" - is critical for accurate modeling of 
spatial correlations on scales > 10% of simulation box size. We provide further illustration of the 
effects of the DC mode on the abundance of halos in small boxes and show that it is straightforward 
to incorporate this mode in cosmological codes that use the "supercomoving" variables. The equations 
governing evolution of dark matter and baryons recast with these variables are particularly simple 
and include the expansion factor, and hence the effect of the DC mode, explicitly only in the Poisson 
equation. 

Subject headings: cosmology: theory - methods: numerical 
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1. INTRODUCTION 

Cosmological simulations have become the main the- 
oretical tool for studying the evolution of cosmic struc- 
tures. Their applications range from modeling the sub- 
parsec scale environments of first stars and supermas- 
sive black holes to the large-scale distribution of matter 
and galaxies. Correspondingly, sizes of simulated regions 
range from hundreds of kiloparsecs to gigaparsecs. 

Many modern simulations use small simulation vol- 
umes, focusing limited computational resources on re- 
solving important small-scale dynamics of galaxies or 
dark matter substructure. In order for such a simulation 
to remain a fair realization of a region of the universe, it 
must properly account for the time evolution of the non- 
vanishing fluctuation of the cosmi c dens i ty at the scale 
of the simulation box. Following iSirkol (|2005[ ) , we call 
this fluctuation the "DC mode,"Q £>dc, to reflect the fact 
that it is constant in space over the simulation volume. 

While any simulation of a finite volume has a non- 
vanishing DC mode, it is not easy to compute for a vol- 
ume of arbitrary geometry. Therefore, in the following 
we assume that the simulation volume is a cubic box 
of size L, and that periodic boundary conditions are im- 
posed along each box dimension. This is indeed the most 
common setup for cosmological simulations. 

Throughout this paper we use a cosmological model 
consistent with the third year WMAP results (fi m 
0.24. ft 100 = 0.73, cr 8 = 0.75, and n s = 0.95), as this 
is the model used in our reference 80/i _1 Mpc simulation 
that we adopt as an approximation to a representative 
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volume of the universe. None of our conclusions, how- 
ever, are dependent on the specific values of cosmological 
parameters. 

2. EFFECTS OF THE DC MODE IN COSMOLOGICAL 
SIMULATIONS 

In the infinite universe the power spectrum of density 
fluctuations, P(k), and their correlation function, £(r), 
are the Fourier transforms of each other. In a cubic pe- 
riodic simulation box of finite size this is no longer true. 
Since the real space quantities are directly related to ob- 
servables, one can argue that it is more important to 
maintai n correct £ (r) than P(k) within the simulation 
volume ()Penlll997l) . In this case the power spectrum of 
the density fluctuations inside the simulation box of size 
L becomes a convolution of the true cosmic power spec- 
trum P{k) and the window function of the simulation 

volume Wl{Ic), 



^grid(fc) = / P(k')W L (k-k')d D k', 



(1) 



where 



D 



and D is the dimension of space. Even if P(0) = 0, 

-Pgrid(O) is not, in general, equal to zero because rms value 
of density fluctuation at the box scale is not zero. Sim- 
ulations aiming to model density fluctuations on scales 
comparable to box size correctly must therefore incorpo- 
rate a non-zero DC mode. 

In order to illustrate the effect of using P gr id instead 
of the true linear P(k) while generating the cosmologi- 
cal initial conditions, we create a large number of real- 
izations of the linear density field for each value of the 
simulation box size L and measure the linear mass vari- 
ance on the scale of 8/i~ 4 Mpc, eg, as the average over 
the whole ensemble^. Figure [T] shows erg as a function of 

8 The specific number of realizations, N, for each value of L is 
determined by the requirement that the average <rg is measured to 
1% precision. 
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Fig. 1. — The linear mass variance in spheres of radius 8h~ 1 Mpc 
(erg) at z = as a function of the simulation box size L, as measured 
from a large number of realizations. Red dotted line shows the 
standard method, where the true linear P(k) is used to generate 
the cosmic density field, while the blue solid line shows the method 
that uses P gr id for setting the density amplitude on the grid. 

simulation box size L at z = for two sets of realizations 
- one that uses P gr id and another one that uses true lin- 
ear P(k). In the former method the linear mass variance 
in spheres is preserved to at least 90% at scales as small 
as half the box size, while in the stan dard method (e.g., 
IBertschingcr 200l t IPrunet et al.ll2008r i the mass variance 
is only accurate in spheres that are less than 10% of the 
box size in radius. The same is true for many other real- 
space clust erin g measures, a s is amply demonstrated by 
IPenl ffT997h and lSirkol d20"05h . 

Using Pg r id to set up initial conditions clearly increases 
the fidelity of a small-box cosmological simulation. How- 
ever, how important is it to have a non-zero DC mode? 
After all, the non-zero value of P gr id(0) only implies a 
non-vanishing rms of a Gaussian-distributed DC mode. 
A single simulation can always impose the constraint of 
Sdc = without violating the statistical properties of 
the initial conditions, although, strictly speaking, such 
initial conditions will not be a true random realization of 
the universe. Hence, if an ensemble of simulations is per- 
formed, then proper sampling of the DC mode is crucial, 
even if the ensemble only contains just two simulations. 

To illustrate this point, we show in Figure [2] the rms 
variation in the number of dark matter halos as a func- 
tion of their maximum circular velocity or the virial 
mass. To compute the rms, we generate two sets of five 
different random realizations of initial conditions, in a 
L = 20/i~ 1 Mpc box. In the first set of initial condi- 
tions the DC mode is forced to be zero, while in the sec- 
ond set the DC mode is computed properly as a Gaus- 
sian distributed random number with the rms value of 
Pgrid (0) / L 3 . Each of the ten realizations of initial condi- 
tions is then evolved to z — with the N- body part of the 
Adaptive Refinement Tree (ART) code (jKravtsov et al.l 
fl997t IKravtsovl [l999h with 128 3 particles and 2 12 dy- 
namic range. 

As a control sample, we use a halo population from a 
single simulation of L — 80/t" 1 Mpc with zero DC mode. 
The simulation followed 512 3 particles with a dynamic 
spatial range of 2 15 . Simulation results of this larger 
box at z = are used to compute the rms number den- 
sity fluctuations of dark matter halos in cubic boxes of 
20/i _1 Mpc on a side. In both simulations halos were 
identified using a variant of the Bound Density Maxima 




Fig. 2. — The rms fluctuations in the numbers of dark matter 
halos as a function of halo maximum circular velocity (top) or 
the virial mass (bottom) in cells of L = 20h~ 1 Mpc. Dotted blue 
and dashed red lines show the average over two ensembles of 5 
simulations each, with the DC mode forced to zero and the properly 
computed DC mode respectively. Simulations followed evolution in 
a box of L = 20h~ 1 Mpc; each of the ten simulations used different 
random realizations of initial conditions. Solid black lines show the 
rms computed from a single simulation of the same cosmology of 
L = 80/i _1 Mpc box size. 
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Fig. 3. — The rms amplitude of the DC mode at z 
function of the simulation box size L. 
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algorithm, as described in iKravtsov et al.l ((20041 ) . 

As Figure [2] illustrates, the ensemble of simulations 
with the DC mode (red dashed lines) properly accounts 
for the fluctuations on the box size scale. Setting the 
DC mode to zero (blue dotted lines) results instead in a 
factor of ~ 2 — 3 underestimate of the true variance. 

3. INCORPORATING THE DC MODE IN A 
COSMOLOGICAL CODE 

The DC mode is substantial for many of the commonly 
used simulation box sizes. For example, Figure [3] shows 
the rms DC mode as a function of the simulation box 
size. The rms DC mode falls below 1% only for box sizes 
L > 500/i _1 Mpc. While 1% may seem like a small num- 
ber, many studies use such boxes to calibrate statistics, 
such as power spectrum and abundance and clustering 
of halos, to accuracy of < 5%. Therefore, effect of the 
DC mode may need to be evaluated even for boxes of 
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hundreds of megaparsecs in size. 

The DC mode thus cannot be "approximately ne- 
glected," at least not in any study aiming at obtaining 
correct statistics on the scale of the simulation box size. 
How easy it is to incorporate the DC mode depends on 
the nature of a cosmological simulation code and could 
be a non-trivial task. 

Cosmological simulations most commonly employ pe- 
riodic boundary conditions and simulation box can be 
considered as a separate universe that expands at a dif- 
ferent rate than the target model universe. Following 
ISirkol (pOOl . we define two expansion factors: the scale 
factor of the simulation box, a^ox, and the true, global 
expansion factor of the universeas, a unl . As mass is con- 
served, the average mass density in the simulation box 
should be equal to the density of a region with the over- 
density odc in the target model universe, 



= -3— [1 + dDc(a U ni)J 



a box °uni 



or 



[l+A DC £>+(auni)] 



1/3 ' 



(2) 
(3) 



where we explicitly spelled out the dependence of the DC 
mode on the universal expansion factor a un i, #Dc(auni) = 
ADC-D+(auni)- Here D + is the linear growth factor of 
density perturbations; it is one (namely, growing) of 
the two independent solutions of the linear perturbation 
equation for d ust-like matter in the Newt onian approxi- 
mation^] (e.g., lBonnorl[l957t lPeebleslfl980T ): 



6 + 2-8 = 4nGp6, 
a 



(4) 



where 6 = (p — p)/p is the matter overdensity with re- 
spect to the mean density of the universe p, a(t) is the 
expansion factor, and derivatives are taken with respect 
to the physical time t. 

We adopt a normalization for D+ such that in a uni- 
verse filled with matter and radiation only (i.e., at suffi- 
ciently early times), 



D + (a)=a + -a eq 



'<•<! 



2VT 



21n(2) 



hi 



VTTx-i 



Vi + x + 1 



. (5) 



9 One may question the validity of the Newtonian approxima- 
tion on large scales in simulations of large box sizes, L = O(c/Ho), 
given that evolution of large-scale modes should be governed by 
fully relativistic equations. However, we note that equation (QJ 
is identical to the corresponding fully relativistic linear perturba- 
tion growth equation in the synchronous gauge or, more generally, 
in the comoving "total matter gauge". The equations governing 
evolution of pert urbations in the synchrono us gauge are presented, 
for example, in Ma & Bcrtschingcr (1995). We can derive equa- 
tion identical to eq. l|4j above by 1) substituting their eq. 21a 
into 21c, 2) considering the case of universe with energy density 
dominated by non- relativistic matter (<5T ° = pS, 5T? = 0), for 
which dh/dr = —26 (their eq. 42), and 3) by making appropriate 
transformations from conformal time r to physical time t using 
dr = dt/a(r). This means that linear evolution of perturbations 
is treated correctly in the Newtonian approximation of cosmolog- 
ical codes at any scale, provided one inte rprets results of simu- 
lations in the appropriate g auge (see, e.g.. lWands &; Slosarl 120091 ; 
IChisari fe Zaldarriaga 20lJ). 



where x = a/a eq and a eq = CIr/0.m is the scale factor of 
matter-radiaton equality. The last two terms can often 
be neglected, the third one falls below 1% for z < 1,100 
a nd the secon d one is below 1% for z < 35. 

ISirkol (|2005D also presents a different form for the a box — 
a U ni relation, which he calls a "Lagrangian viewpoint" as 
opposite to the "Elurian viewpoint" of equation (J3|) , but 
in essence is just a first order expansion of Equation ([3]) 
in A DC , 



abc 



1 - ^A DC £> + (a uni ) 



(6) 



When we use this to compute the rms variation in the 
number of dark matter halos shown in Figure [21 we find a 
result which is virtually indistinguishable from the "Eule- 
rian viewpoint" shown with red dashed lines. We choose 
to use the "Eulerian viewpoint" as our primary method 
for accounting for the DC mode, because it explicitly 
conserves mass to all orders in the perturbation theory. 

The parameter Adc is constant in a given simulation 
and describes the amplitude of the DC mode in a given 
realization of initial conditions. In principle, equation ([3]) 
is also valid in the non-linear regime, if we treat D+ as 
a specific non-linear growth mode in a given simulation 
box that properly accounts for the coupling of the DC 
mode to other modes (including modes with wavelengths 
larger than L). In practice, however, there exists no way 
to compute the non-linear growth rate besides the nu- 
merical simulation itself, so hereafter we assume that the 
DC mode remains with sufficient precision in the linear 
regime throughout the whole time-span of the numerical 
simulation. Of course, the term "sufficient precision" de- 
pends on the required precision of the simulation results. 
For example, for our ensembles of 20/i _1 Mpc boxes the 
rms value of odc(1) is almost 0.7; never-the-less, the 
linear approximation for the DC mode is sufficient to 
achieve about 20% agreement between an ensemble of 
20ra _1 Mpc simulations and a single 80/i _1 Mpc run in 
Figure El 

Historically, the DC mode was sometimes incorporated 
into a simulation by appropriately rescaling the values of 
cosmological parameters. An examp l e of s uch rescaling 
is also given in the appendix of lSirkol (|2005l ). Effectively, 
in a spatially flat cosmology such a rescaling introduces 
non-zero spatial curvature (or additional curvature if the 
universe is assumed to be spatially non-flat) with the 
extra curvature parameter AQk ~ — (5/3)fijif Adc- 

However, such approach does not work in more general 
cosmological models. For example, one cannot consider 
simulation with a DC mode as a separate universe in 
models with non-negligible amount of relativistic matter, 
general dark energy component, decaying dark matter, or 
modified gravity. 

We, therefore, advocate an alternative approach that 
can be used to include the DC mode in simulation codes 
exactly. This can be done by following both abox and a un i 
as a function of cosmic time and incorporating this differ- 
ence in simulation equations explicitly. A disadvantage 
of this approach is that some of the equations become 
unnecessarily complicated. For example, the evolution 
equation for the peculiar velocity v pec of a dark matter 
particle in the expanding coordinate system commonly 
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used in cosmological codes is 



dv 



pec 



dt 



where dot symbolizes the time derivative, <f> is the pe- 
culiar gravitational potential, and spatial derivatives are 
taken with respect to the comoving coordinates x, which 
is emphasized by the subscript x under the gradient sym- 
bol. Computing ctbox requires a differentiation of equa- 
tion ([3]) and leads to uncouth and inelegant (but fully 
tractable) formulae. 

A much more elegant way to incorporate the 
DC mode in a cosmological code is via the so- 
called supercomovi n g var iableS^\ first introduced by 
iDoroshkevich et al~l (11980T). and then us ed in several 
nume rica l works ( Shandarinl Il980l: IShapiro et al. 
1981 IShapiro fc Struck-Marcelll Il985t IGnedin 



19951: lYepes et all 119971: IMartel fc Shapiro! 119981: 
Kravtsov et all l2002HTevssierl |200^T Specifically, 
the proper coordinates r, velocity u, density p, and 
gravitational potential $ are replaced with comoving 
coordinates x, supercomoving peculiar velocities v (to be 
distinguished from the normal peculiar velocities v pec ), 
comoving density g, and peculiar supercomoving grav- 
itational potential <p = a 2 4> according to the following 
transformations: 



= ax, 



1 



u = aHx H — v, 
a 

Q 

P=-T, 

a 6 

ad o 1 
$ = — — x z + -y(p. 



(7) 
(8) 
(9) 
(10) 



In addition, the cosmic time t is replaced with the super- 
comoving time, r, defined such that 



dt 

dr = —x. 

er 



(11) 



In variables t, x, v, g, (p the equations of motion of dark 
matter particles assume the form identical to that in the 
proper, non-expanding reference frame, without any ex- 
plicit cosmological terms, 

dx 
dv 

Similarly, the Euler equations for monatomic gas (poly- 
tropic index 7 = 5/3) include no explicit cosmological 
terms either. The only equation, in which the cosmic 
scale factor a appears explicitly, is the Poisson equation, 



Vlip = 4irGa(g - g) = AirGagS, 



(12) 



where g is the mean comoving matter density of the uni- 
verse. Note that g = const in cosmological models in 
which the dark matter does not decay into radiation. 

Incorporating the DC mode into a simulation code us- 
ing the supercomoving variables is straightforward: all 

10 T he term "supercomoving" was coined by[Martcl & Shapiro] 
CE99|). 



that is required is to identify the scale factor a in equa- 
tion (|12p with the box scale factor abox, 



\?lip = 4nGa ho ^gS, 



(13) 



and compute both abox and a un ; as a function of the 
supercomoving time r. 

For simulations including gas and other physics, abox 
should be used anywhere a code would use the scale fac- 
tor internally. For example, in unit conversions, in cool- 
ing functions and star formation and feedback recipes, 
for checking energy conservation, in the radiative trans- 
fer solver, etc. For some physical processes, for example 
redshift-dependent cosmic backgrounds (radiation, cos- 
mic rays, etc), the proper choice of scaling is model 
dependent. If a hydrodynamic simulation uses cooling 
rates that account for the cosmic ultraviolet background 
in a t abulated form (e.g.. lKravtsovl[2003l : iWiersma et al.1 
2009), the table should typically be evaluated at a un j. 
If the sources of the radiation are inhomogeneous on the 
scale of the simulation, however, evaluating at abox would 
more accurately capture the large-scale density depen- 
dence. The only remaining use for a un i is in identifying 
the snapshot epoch and the final epoch (if a simulation 
is evolved to z = 0, the stopping criterion is a unl = 1, 
not abox = 1). 

4. CONCLUSIONS 

In this short note we illustrated a known but not widely 
appreciated fact that taking the DC mode into account 
in small box cosmological simulations is a requirement 
for the simulation to serve as a representative volume of 
the universe and to model density fluctuations correctly. 
The DC mode must be included in simulations if more 
than a single realization of the initial conditions is used. 

While approximate methods for including the DC 
mode have been proposed, here we advocate the use of 
the supercomoving variables in cosmological simulations 
with which the effect of the DC mode is limited to a sim- 
ple multiplicative term in the Poisson equation (|13|) and 
can be incorporated simply and exactly. 

An efficient and reliable cosmology module that com- 
putes several important quantities (abox, a un - u superco- 
moving time r, the linear growth rate D +1 etc.) via 
lookup tables in the linear regime is available from the 
authors upon request. 
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